subroutine sort_and_discard()

  use logs
  use indexx
  use particle_data
  implicit none
  integer*4 :: ipart,itype,i0,i1

  write(logunit,*)' -performing global sort...'
  if(NumPart_Total_AllSnaps(0) > 0) call indexxx(NumPart_Total_AllSnaps(0),ParticleIDs%Type0,IndexIDs%Type0)
  if(NumPart_Total_AllSnaps(1) > 0) call indexxx(NumPart_Total_AllSnaps(1),ParticleIDs%Type1,IndexIDs%Type1)
  if(NumPart_Total_AllSnaps(2) > 0) call indexxx(NumPart_Total_AllSnaps(2),ParticleIDs%Type2,IndexIDs%Type2)
  if(NumPart_Total_AllSnaps(3) > 0) call indexxx(NumPart_Total_AllSnaps(3),ParticleIDs%Type3,IndexIDs%Type3)
  if(NumPart_Total_AllSnaps(4) > 0) call indexxx(NumPart_Total_AllSnaps(4),ParticleIDs%Type4,IndexIDs%Type4)
  if(NumPart_Total_AllSnaps(5) > 0) call indexxx(NumPart_Total_AllSnaps(5),ParticleIDs%Type5,IndexIDs%Type5)

  ! Shift index lists into range [0:n-1]
  if(NumPart_Total_AllSnaps(0) > 0) IndexIDs%Type0(:) = IndexIDs%Type0(:) - 1
  if(NumPart_Total_AllSnaps(1) > 0) IndexIDs%Type1(:) = IndexIDs%Type1(:) - 1
  if(NumPart_Total_AllSnaps(2) > 0) IndexIDs%Type2(:) = IndexIDs%Type2(:) - 1
  if(NumPart_Total_AllSnaps(3) > 0) IndexIDs%Type3(:) = IndexIDs%Type3(:) - 1
  if(NumPart_Total_AllSnaps(4) > 0) IndexIDs%Type4(:) = IndexIDs%Type4(:) - 1
  if(NumPart_Total_AllSnaps(5) > 0) IndexIDs%Type5(:) = IndexIDs%Type5(:) - 1

  ! Fix properties of first particle for each type manually
  do itype=0,5
     if(NumPart_Total_AllSnaps(itype) > 0) NumPart_Unique(itype) = 1
  enddo
  if(NumPart_Total_AllSnaps(0) > 0) UniqueIDs%Type0(0)     = ParticleIDs%Type0(IndexIDs%Type0(0))
  if(NumPart_Total_AllSnaps(1) > 0) UniqueIDs%Type1(0)     = ParticleIDs%Type1(IndexIDs%Type1(0))
  if(NumPart_Total_AllSnaps(2) > 0) UniqueIDs%Type2(0)     = ParticleIDs%Type2(IndexIDs%Type2(0))
  if(NumPart_Total_AllSnaps(3) > 0) UniqueIDs%Type3(0)     = ParticleIDs%Type3(IndexIDs%Type3(0))
  if(NumPart_Total_AllSnaps(4) > 0) UniqueIDs%Type4(0)     = ParticleIDs%Type4(IndexIDs%Type4(0))
  if(NumPart_Total_AllSnaps(5) > 0) UniqueIDs%Type5(0)     = ParticleIDs%Type5(IndexIDs%Type5(0))
  if(NumPart_Total_AllSnaps(0) > 0) UniqueSnapIDs%Type0(0) = SnapIDs%Type0(IndexIDs%Type0(0))
  if(NumPart_Total_AllSnaps(1) > 0) UniqueSnapIDs%Type1(0) = SnapIDs%Type1(IndexIDs%Type1(0))
  if(NumPart_Total_AllSnaps(2) > 0) UniqueSnapIDs%Type2(0) = SnapIDs%Type2(IndexIDs%Type2(0))
  if(NumPart_Total_AllSnaps(3) > 0) UniqueSnapIDs%Type3(0) = SnapIDs%Type3(IndexIDs%Type3(0))
  if(NumPart_Total_AllSnaps(4) > 0) UniqueSnapIDs%Type4(0) = SnapIDs%Type4(IndexIDs%Type4(0))
  if(NumPart_Total_AllSnaps(5) > 0) UniqueSnapIDs%Type5(0) = SnapIDs%Type5(IndexIDs%Type5(0))

  write(logunit,*)' -discarding duplicates...'
  itype = 0
  do ipart=1,NumPart_Total_AllSnaps(itype)-1
     i1 = IndexIDs%Type0(ipart)
     i0 = IndexIDs%Type0(ipart-1)
     if(ParticleIDs%Type0(i1) /= ParticleIDs%Type0(i0))then 
        ! Unique particle, so add to list 
        UniqueIDs%Type0(NumPart_Unique(itype))     = ParticleIDs%Type0(i1)
        UniqueSnapIDs%Type0(NumPart_Unique(itype)) = SnapIDs%Type0(i1)
        NumPart_Unique(itype)                      = NumPart_Unique(itype) + 1
     else
        ! Duplicate particle, so just need to check whether 
        ! this copy came from a later snapshot
        if(SnapIDs%Type0(i1) > UniqueSnapIDs%Type0(NumPart_Unique(itype)-1)) then
           UniqueSnapIDs%Type0(NumPart_Unique(itype)-1) = SnapIDs%Type0(i1)
        endif
     endif
  enddo

  itype = 1
  do ipart=1,NumPart_Total_AllSnaps(itype)-1
     i1 = IndexIDs%Type1(ipart)
     i0 = IndexIDs%Type1(ipart-1)
     if(ParticleIDs%Type1(i1) /= ParticleIDs%Type1(i0))then 
        ! Unique particle, so add to list 
        UniqueIDs%Type1(NumPart_Unique(itype))     = ParticleIDs%Type1(i1)
        UniqueSnapIDs%Type1(NumPart_Unique(itype)) = SnapIDs%Type1(i1)
        NumPart_Unique(itype)                      = NumPart_Unique(itype) + 1
     else
        ! Duplicate particle, so just need to check whether 
        ! this copy came from a later snapshot
        if(SnapIDs%Type1(i1) > UniqueSnapIDs%Type1(NumPart_Unique(itype)-1)) then
           UniqueSnapIDs%Type1(NumPart_Unique(itype)-1) = SnapIDs%Type1(i1)
        endif
     endif
  enddo

  itype = 2
  do ipart=1,NumPart_Total_AllSnaps(itype)-1
     i1 = IndexIDs%Type2(ipart)
     i0 = IndexIDs%Type2(ipart-1)
     if(ParticleIDs%Type2(i1) /= ParticleIDs%Type2(i0))then 
        ! Unique particle, so add to list 
        UniqueIDs%Type2(NumPart_Unique(itype))     = ParticleIDs%Type2(i1)
        UniqueSnapIDs%Type2(NumPart_Unique(itype)) = SnapIDs%Type2(i1)
        NumPart_Unique(itype)                      = NumPart_Unique(itype) + 1
     else
        ! Duplicate particle, so just need to check whether 
        ! this copy came from a later snapshot
        if(SnapIDs%Type2(i1) > UniqueSnapIDs%Type2(NumPart_Unique(itype)-1)) then
           UniqueSnapIDs%Type2(NumPart_Unique(itype)-1) = SnapIDs%Type2(i1)
        endif
     endif
  enddo

  itype = 3
  do ipart=1,NumPart_Total_AllSnaps(itype)-1
     i1 = IndexIDs%Type3(ipart)
     i0 = IndexIDs%Type3(ipart-1)
     if(ParticleIDs%Type3(i1) /= ParticleIDs%Type3(i0))then 
        ! Unique particle, so add to list 
        UniqueIDs%Type3(NumPart_Unique(itype))     = ParticleIDs%Type3(i1)
        UniqueSnapIDs%Type3(NumPart_Unique(itype)) = SnapIDs%Type3(i1)
        NumPart_Unique(itype)                      = NumPart_Unique(itype) + 1
     else
        ! Duplicate particle, so just need to check whether 
        ! this copy came from a later snapshot
        if(SnapIDs%Type3(i1) > UniqueSnapIDs%Type3(NumPart_Unique(itype)-1)) then
           UniqueSnapIDs%Type3(NumPart_Unique(itype)-1) = SnapIDs%Type3(i1)
        endif
     endif
  enddo

  itype = 4
  do ipart=1,NumPart_Total_AllSnaps(itype)-1
     i1 = IndexIDs%Type4(ipart)
     i0 = IndexIDs%Type4(ipart-1)
     if(ParticleIDs%Type4(i1) /= ParticleIDs%Type4(i0))then 
        ! Unique particle, so add to list 
        UniqueIDs%Type4(NumPart_Unique(itype))     = ParticleIDs%Type4(i1)
        UniqueSnapIDs%Type4(NumPart_Unique(itype)) = SnapIDs%Type4(i1)
        NumPart_Unique(itype)                      = NumPart_Unique(itype) + 1
     else
        ! Duplicate particle, so just need to check whether 
        ! this copy came from a later snapshot
        if(SnapIDs%Type4(i1) > UniqueSnapIDs%Type4(NumPart_Unique(itype)-1)) then
           UniqueSnapIDs%Type4(NumPart_Unique(itype)-1) = SnapIDs%Type4(i1)
        endif
     endif
  enddo

  itype = 5
  do ipart=1,NumPart_Total_AllSnaps(itype)-1
     i1 = IndexIDs%Type5(ipart)
     i0 = IndexIDs%Type5(ipart-1)
     if(ParticleIDs%Type5(i1) /= ParticleIDs%Type5(i0))then 
        ! Unique particle, so add to list 
        UniqueIDs%Type5(NumPart_Unique(itype))     = ParticleIDs%Type5(i1)
        UniqueSnapIDs%Type5(NumPart_Unique(itype)) = SnapIDs%Type5(i1)
        NumPart_Unique(itype)                      = NumPart_Unique(itype) + 1
     else
        ! Duplicate particle, so just need to check whether 
        ! this copy came from a later snapshot
        if(SnapIDs%Type5(i1) > UniqueSnapIDs%Type5(NumPart_Unique(itype)-1)) then
           UniqueSnapIDs%Type5(NumPart_Unique(itype)-1) = SnapIDs%Type5(i1)
        endif
     endif
  enddo

  ! Free up some memory
  deallocate(ParticleIDs%Type0)
  deallocate(ParticleIDs%Type1)
  deallocate(ParticleIDs%Type2)
  deallocate(ParticleIDs%Type3)
  deallocate(ParticleIDs%Type4)
  deallocate(ParticleIDs%Type5)
  !
  deallocate(IndexIDs%Type0)
  deallocate(IndexIDs%Type1)
  deallocate(IndexIDs%Type2)
  deallocate(IndexIDs%Type3)
  deallocate(IndexIDs%Type4)
  deallocate(IndexIDs%Type5)
  !
  deallocate(SnapIDs%Type0)
  deallocate(SnapIDs%Type1)
  deallocate(SnapIDs%Type2)
  deallocate(SnapIDs%Type3)
  deallocate(SnapIDs%Type4)
  deallocate(SnapIDs%Type5)




end subroutine sort_and_discard



subroutine match_to_ics()

 use logs
 use indexx
 use particle_data
 use ic_data

 implicit none
 integer*4 :: ipart,itype,jpart
 integer*4 :: id_seek,id_current
 
 write(logunit,*)' -sorting ParticleIDs from the ICs'
 call flush(logunit)

 ! Sort the ics ids - note that we only need to sort those types we'll be matching against
 if(NumPart_ICs(0) > 0) then 
    write(logunit,*)' -Type0'
    call indexxx(NumPart_ICs(0),IC_ParticleIDs%Type0,IC_IndexIDs%Type0) 
    IC_IndexIDs%Type0(:) = IC_IndexIDs%Type0(:) -1
 endif
 if(NumPart_ICs(1) > 0) then  
    write(logunit,*)' -Type1'
    call indexxx(NumPart_ICs(1),IC_ParticleIDs%Type1,IC_IndexIDs%Type1) 
    IC_IndexIDs%Type1(:) = IC_IndexIDs%Type1(:) -1
 endif
 if(NumPart_ICs(2) > 0) then 
    write(logunit,*)' -Type2'
    call indexxx(NumPart_ICs(2),IC_ParticleIDs%Type2,IC_IndexIDs%Type2)
    IC_IndexIDs%Type2(:) = IC_IndexIDs%Type2(:) -1
 endif
 if(NumPart_ICs(3) > 0) then 
    write(logunit,*)' -Type3'
    call indexxx(NumPart_ICs(3),IC_ParticleIDs%Type3,IC_IndexIDs%Type3)
    IC_IndexIDs%Type3(:) = IC_IndexIDs%Type3(:) -1
 endif
#ifdef TYPE4_IS_CDM
 if(NumPart_ICs(4) > 0) then 
    write(logunit,*)' -Type4'
    call indexxx(NumPart_ICs(4),IC_ParticleIDs%Type4,IC_IndexIDs%Type4)
    IC_IndexIDs%Type4(:) = IC_IndexIDs%Type4(:) -1
 endif
#endif
 if(NumPart_ICs(5) > 0) then 
    write(logunit,*)' -Type5'
    call indexxx(NumPart_ICs(5),IC_ParticleIDs%Type5,IC_IndexIDs%Type5)
    IC_IndexIDs%Type5(:) = IC_IndexIDs%Type5(:) -1
 endif

 ! For each unique selected ID, match to the ics id 
 write(logunit,*)' -matching selected particles to ICs'
 call flush(logunit)

 ! The procedure here depends on the type of ICs
 ! Gas is often added by the code, and so is absent from the 
 ! ICs - in this case, we match to the partner CDM particle
 ! Also, selected stars will be gas by the time of the ICs, so 
 ! if gas is present, we match to this, if not, we match to CDM.

 write(logunit,*)' -Type0'
 call flush(logunit)
 if(NumPart_ICs(0) > 0)then ! Regular behaviour

    write(logunit,*)'  -matching Type0 to IC Type0'

    jpart = 0
    do ipart=0,NumPart_Unique(0)-1     
       id_seek    = UniqueIDs%Type0(ipart)
       id_current = IC_ParticleIDs%Type0(IC_IndexIDs%Type0(jpart))
       do while(id_current < id_seek)
          if(jpart == NumPart_ICs(0)) call abort(" -cannot match selected Type0 id to ICs ID")
          jpart = jpart + 1
          id_current = IC_ParticleIDs%Type0(IC_IndexIDs%Type0(jpart))
       end do
       UniqueCoordinates%Type0(0:2,ipart) = IC_Coordinates%Type0(0:2,IC_IndexIDs%Type0(jpart))
       UniqueVelocity%Type0(0:2,ipart)    = IC_Velocity%Type0(0:2,IC_IndexIDs%Type0(jpart))
    enddo
    
 else !Match associated CDM particle to Type1
    
    write(logunit,*)'  -matching Type0 to IC associated Type1'

    jpart = 0
    do ipart=0,NumPart_Unique(0)-1     
       id_seek    = UniqueIDs%Type0(ipart) - 1000000000 ! Converts gas ID to associated CDM ID
       id_current = IC_ParticleIDs%Type1(IC_IndexIDs%Type1(jpart))
       do while(id_current < id_seek)
          if(jpart == NumPart_ICs(1)) call abort(" -cannot match selected Type0 id to associated CDM ICs ID")
          jpart = jpart + 1
          id_current = IC_ParticleIDs%Type1(IC_IndexIDs%Type1(jpart))
       end do
       UniqueCoordinates%Type0(0:2,ipart) = IC_Coordinates%Type1(0:2,IC_IndexIDs%Type1(jpart))
       UniqueVelocity%Type0(0:2,ipart)    = IC_Velocity%Type1(0:2,IC_IndexIDs%Type1(jpart))
    enddo
    
 endif

 ! Types 1-3, 5 are relatively simple....
 write(logunit,*)'  -matching Type1 to IC Type1'
 itype = 1 
 jpart = 0
 do ipart=0,NumPart_Unique(itype)-1     
    id_seek    = UniqueIDs%Type1(ipart)
    id_current = IC_ParticleIDs%Type1(IC_IndexIDs%Type1(jpart))
    do while(id_current < id_seek)
       if(jpart == NumPart_ICs(itype)) call abort(" -cannot match selected Type1 id to ICs ID")
       jpart = jpart + 1
       id_current = IC_ParticleIDs%Type1(IC_IndexIDs%Type1(jpart))
    end do
    UniqueCoordinates%Type1(0:2,ipart) = IC_Coordinates%Type1(0:2,IC_IndexIDs%Type1(jpart))
    UniqueVelocity%Type1(0:2,ipart)    = IC_Velocity%Type1(0:2,IC_IndexIDs%Type1(jpart))
 enddo

 write(logunit,*)'  -matching Type2 to IC Type2'
 itype = 2 
 jpart = 0
 do ipart=0,NumPart_Unique(itype)-1     
    id_seek    = UniqueIDs%Type2(ipart)
    id_current = IC_ParticleIDs%Type2(IC_IndexIDs%Type1(jpart))
    do while(id_current < id_seek)
       if(jpart == NumPart_ICs(itype)) call abort(" -cannot match selected Type2 id to ICs ID")
       jpart = jpart + 1
       id_current = IC_ParticleIDs%Type2(IC_IndexIDs%Type2(jpart))
    end do
    UniqueCoordinates%Type2(0:2,ipart) = IC_Coordinates%Type2(0:2,IC_IndexIDs%Type2(jpart))
    UniqueVelocity%Type2(0:2,ipart)    = IC_Velocity%Type2(0:2,IC_IndexIDs%Type2(jpart))
 enddo

 write(logunit,*)'  -matching Type3 to IC Type3'
 itype = 3 
 jpart = 0
 do ipart=0,NumPart_Unique(itype)-1     
    id_seek    = UniqueIDs%Type3(ipart)
    id_current = IC_ParticleIDs%Type3(IC_IndexIDs%Type3(jpart))
    do while(id_current < id_seek)
       if(jpart == NumPart_ICs(itype)) call abort(" -cannot match selected Type3 id to ICs ID")
       jpart = jpart + 1
       id_current = IC_ParticleIDs%Type3(IC_IndexIDs%Type3(jpart))
    end do
    UniqueCoordinates%Type3(0:2,ipart) = IC_Coordinates%Type3(0:2,IC_IndexIDs%Type3(jpart))
    UniqueVelocity%Type3(0:2,ipart)    = IC_Velocity%Type3(0:2,IC_IndexIDs%Type3(jpart))
 enddo

#ifdef TYPE4_IS_CDM
 ! Simple case, type4 is just another collissionless species
 write(logunit,*)'  -matching Type4 to IC Type4'
 itype = 4 
 jpart = 0
 do ipart=0,NumPart_Unique(itype)-1     
    id_seek    = UniqueIDs%Type4(ipart)
    id_current = IC_ParticleIDs%Type4(IC_IndexIDs%Type4(jpart))
    do while(id_current < id_seek)
       if(jpart == NumPart_ICs(itype)) call abort(" -cannot match selected Type4 id to ICs ID")
       jpart = jpart + 1
       id_current = IC_ParticleIDs%Type4(IC_IndexIDs%Type4(jpart))
    end do
    UniqueCoordinates%Type4(0:2,ipart) = IC_Coordinates%Type4(0:2,IC_IndexIDs%Type4(jpart))
    UniqueVelocity%Type4(0:2,ipart)    = IC_Velocity%Type4(0:2,IC_IndexIDs%Type4(jpart))
 enddo

#else
 ! Type4 are stars, so must be matched to gas, if it is present, or the 
 ! CDM associated with the parent gas particle from whence each star came. Easy.

 if(NumPart_ICs(0) > 0)then ! Associate stars to parent gas particle

    write(logunit,*)'  -matching Type4 to IC Type0'
    jpart = 0
    do ipart=0,NumPart_Unique(4)-1     
       id_seek    = UniqueIDs%Type4(ipart)
       id_current = IC_ParticleIDs%Type0(IC_IndexIDs%Type0(jpart))
       do while(id_current < id_seek)
          if(jpart == NumPart_ICs(0)) call abort(" -cannot match selected Type4 id to parent gas ICs ID")
          jpart = jpart + 1
          id_current = IC_ParticleIDs%Type0(IC_IndexIDs%Type0(jpart))
       end do
       UniqueCoordinates%Type4(0:2,ipart) = IC_Coordinates%Type0(0:2,IC_IndexIDs%Type0(jpart))
       UniqueVelocity%Type4(0:2,ipart)    = IC_Velocity%Type0(0:2,IC_IndexIDs%Type0(jpart))
    enddo
    
 else ! Match parent gas particle to associated CDM particle in Type1 
    
    write(logunit,*)'  -matching Type4 to IC associated Type1'
    jpart = 0
    do ipart=0,NumPart_Unique(4)-1     

! This bit depends upon the assumptions in the ICs - sometimes gas is given its 
! associated CDM particle's ID + 10^9, otherwise the gas IDs just run 1:NumPart_Total(0)
#ifdef GASIDS_ADD_BILLION
       id_seek    = UniqueIDs%Type4(ipart) - 1000000000     ! Converts stellar ID to associated CDM ID of parent gas particle
#else
       id_seek    = UniqueIDs%Type4(ipart) + NumPart_ICs(0) ! Converts stellar ID to associated CDM ID of parent gas particle
#endif
       id_current = IC_ParticleIDs%Type1(IC_IndexIDs%Type1(jpart))
       do while(id_current < id_seek)
          if(jpart == NumPart_ICs(1)) call abort(" -cannot match selected Type4 id to parent gas particle's associated CDM ICs ID")
          jpart = jpart + 1
          id_current = IC_ParticleIDs%Type1(IC_IndexIDs%Type1(jpart))
       end do
       UniqueCoordinates%Type4(0:2,ipart) = IC_Coordinates%Type1(0:2,IC_IndexIDs%Type1(jpart))
       UniqueVelocity%Type4(0:2,ipart)    = IC_Velocity%Type1(0:2,IC_IndexIDs%Type1(jpart))
    enddo
    
 endif

#endif

 write(logunit,*)'  -matching Type5 to IC Type5'
 ! This is simple case again
 itype = 5 
 jpart = 0
 do ipart=0,NumPart_Unique(itype)-1     
    id_seek    = UniqueIDs%Type5(ipart)
    id_current = IC_ParticleIDs%Type5(IC_IndexIDs%Type5(jpart))
    do while(id_current < id_seek)
       if(jpart == NumPart_ICs(itype)) call abort(" -cannot match selected Type5 id to ICs ID")
       jpart = jpart + 1
       id_current = IC_ParticleIDs%Type5(IC_IndexIDs%Type5(jpart))
    end do
    UniqueCoordinates%Type5(0:2,ipart) = IC_Coordinates%Type5(0:2,IC_IndexIDs%Type5(jpart))
    UniqueVelocity%Type5(0:2,ipart)    = IC_Velocity%Type5(0:2,IC_IndexIDs%Type5(jpart))
 enddo

 write(logunit,*)'  -done'
 call flush(logunit)

end subroutine match_to_ics


subroutine reverse_zeldovich
  use logs
  use header_vars
  use particle_data
  implicit none
  integer*4 :: ipart
  real*8    :: HubTerm,prefac

  ! Trace particles back to a -> 0 by reversing the Zel'dovich
  ! approximation
  write(logunit,*)" -applying reverse Zel'dovich transformation"
  write(logunit,*)"   z      = ",binary_ics_header%Redshift

  HubTerm       = (binary_ics_header%Omega0/(binary_ics_header%Time**3.0d0))+binary_ics_header%OmegaLambda
  HubbleParam_a = 1.0d2*sqrt(HubTerm)
  prefac        = real(1.0d0 / (sqrt(binary_ics_header%Time)*HubbleParam_a))

  write(logunit,*)"   H(z)   = ",HubbleParam_a," km/s/Mpc"
  write(logunit,*)"   prefac = ",prefac

  do ipart=0,NumPart_Unique(0)-1
     ForceFreeCoordinates%Type0(0:2,ipart) = UniqueCoordinates%Type0(0:2,ipart) - prefac*UniqueVelocity%Type0(0:2,ipart)
  enddo

  do ipart=0,NumPart_Unique(1)-1
     ForceFreeCoordinates%Type1(0:2,ipart) = UniqueCoordinates%Type1(0:2,ipart) - prefac*UniqueVelocity%Type1(0:2,ipart)
  enddo

  do ipart=0,NumPart_Unique(2)-1
     ForceFreeCoordinates%Type2(0:2,ipart) = UniqueCoordinates%Type2(0:2,ipart) - prefac*UniqueVelocity%Type2(0:2,ipart)
  enddo

  do ipart=0,NumPart_Unique(3)-1
     ForceFreeCoordinates%Type3(0:2,ipart) = UniqueCoordinates%Type3(0:2,ipart) - prefac*UniqueVelocity%Type3(0:2,ipart)
  enddo

  do ipart=0,NumPart_Unique(4)-1
     ForceFreeCoordinates%Type4(0:2,ipart) = UniqueCoordinates%Type4(0:2,ipart) - prefac*UniqueVelocity%Type4(0:2,ipart)
  enddo

  do ipart=0,NumPart_Unique(5)-1
     ForceFreeCoordinates%Type5(0:2,ipart) = UniqueCoordinates%Type5(0:2,ipart) - prefac*UniqueVelocity%Type5(0:2,ipart)
  enddo


end subroutine reverse_zeldovich
